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ABSTRACT 

Interstellar dust grains are responsible for modifying the spectral energy distribution 
(SED) of galaxies, both absorbing starlight at UV and optical wavelengths and con- 
verting this energy into thermal emission in the infrared. The detailed description 
of these phenomena is of fundamental importance in order to compare the predic- 
tions of theoretical models of galaxy formation and evolution with the most recent 
observations in the infrared region. In this paper we compare the results of GRASIL, 
a code explicitly solving for the equation of radiative transfer in a dusty medium, 
with the predictions of a variety of IR template libraries, both analytically and ob- 
servationally determined. We employ star formation history samples extracted from 
the semi-analytical galaxy formation model MORGANA to create libraries of synthetic 
SEDs from the near- to the far-infrared. We consider model predictions at different 
redshift ranges to explore any possible influence in the shape and normalisation of the 
SEDs due to the expected evolution of the galaxy properties. We compute the total 
absorbed starlight predicted by GRASIL at optical wavelengths to statistically compare 
the synthetic SEDs with the selected IR templates. We show that synthetic SEDs at 
a given total infrared luminosity are predicted to be systematically different at dif- 
ferent redshift and for different properties of the underlying model galaxy. However, 
we determine spectral regions where the agreement between the results of radiative 
transfer and IR templates is good in a statistical sense (i.e. in terms of the luminosity 
functions). Moreover, we highlight some potentially relevant discrepancies between the 
different approaches, both in the region dominated by PAH emission and at sub-mm 
wavelengths. These results determine potentially critical issues in the infrared lumi- 
nosity functions as predicted by semi-analytical models coupled with different IR flux 
estimators. 
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1 INTRODUCTION 

The great advance in our understanding of the spectral en- 
ergy distributions (SEDs) of galaxies has clearly demon- 
strated the ubiquitous presence of dust in galaxies. Dust 
plays a major role in galaxy evolution: it contributes to set- 
ting the physical conditions for star formation in galaxies, 
both seeding the formation of molecules and shielding the 
molecular gas from ionising radiation. It is also important in 
the processes of cooling and condensation of giant molecular 
clouds (MCs) to the densities required for the onset of star 



formation fsee lDorschner fc Henninelll995l . for a review). At 
the same time, the presence of dust has a strong impact on 
the intrinsic spectral energy distributions (SED) of galaxies, 
since dust grains efficiently absorb and scatter radiation at 
short wavelengths (A < 1pm). The absorbed energy is then 
thermally re-emitted at longer wavelengths in the infra-red 
(IR) region (A > 1/zm). 

Dust-related processes of attenuation and re-emission 
are key mechanisms in understanding and interpreting galac- 
tic SEDs. Estimates based on the results of the IRAS 
satellite demonstrated that the IR region is responsible for 
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at least ~ 30% of the bolom etric luminosity of nearby 
galaxies llPopescu fc Tuffdl2003 ). Early IR surveys (see, e.g. 
ISanders fc Mirabell 1 19961 7 discovered a population of heav- 
ily obscured luminous and ultra-luminous infrared galax- 
ies (LIRGs, with IR luminosities Lir~ 10 11 — 10 12 L@, and 
ULIRGs, Lm> 10 12 ). Moreover, the study of the cosmic in- 
frared background (see lHauser fe Dwekll200ll . for a review) 
shows that the global energy emitted in the IR is com- 
parable to the direct starlight emission, detectable mainly 
at optical wavelengths, emerging from galaxies. Moreover, 
galaxy number counts in the mid-IR are an order of mag- 
nitude higher than thos e predicted by non - evolving lo- 
cal luminosity functions |Pozzetti et al.l 1 19981; Elbaz et al 
1999 1; lLagache et alJll999h. Several groups dCharv fc Elbaz 



200 J : iDole et al l l200ll ; lLagache et aD l2003h have shown 



that all these pieces of evidence favour strong evolution of 
the bright IR sources. This has been confirmed with di- 
rect o bservations of high-redshift sou rces using facilities like 
ISO l|Dole et al l l200ll ; felbaz et al.l 
2002i) . the Spitzer Sr >ace Telescope 



2002; Gruppioni et al. 


|Le Floc'h et al. 


2005; 


(iChaoman et al. 


2003, 



Babbedee et aD 120061 ) and SCUBA , 
20051 ). This strong evolution of luminous IR sources can also 



account for the peak in the cosmic infrared background. 

Dust reprocessed radiation therefore represents a fun- 
damental aspect of galactic emission: in particular a sig- 
nificant fraction of star formation activity is expected to be 
heavily extinguished and detectable only in the IR, since the 
molecular clouds are both the sites for star formation and 
among the most dusty environments. Most interestingly, this 
contribution was probably more important at higher red- 
shift, where star formation activity and dust reprocessing 
are both stronger. The detailed modeling of dust properties 
is then a key issue in order to understand galaxy evolu- 
tion at different cosmic epochs. Our understanding of these 
processes is complicated by the strong dependence of galac- 
tic SEDs on the relative geometry be tween dust and stars 
|Granato et al.ll2000l ; lTuffs et alj2004 l In fact, the youngest 
stars, which dominate the UV luminosity, are expected to be 
heavily extinguished by their optically thick parent Molec- 
ular Clouds (MC). Older stars survive the disruption of the 
MCs and in fact, an age-dependent extinction has been pos- 
tulate d in the cirrus component for the inte rmediate-age 
stars (jPopescu et al.|[200ol ; IPanuzzo et alj|2007h . in order to 
reproduce the observed properties of disc-dominated galax- 
ies. The emerging SED depends critically on the tempera- 
ture distribution in the different environments (cirrus and 
MCs), as different dust species (with different composition 
and size) respond differently to the radiation field. 

Therefore, a detailed treatment of the effects of dust 
grains in different spectral regions (both attenuation in 
the optical and re-emission in the IR) is of fundamen- 
tal importance in order to infer the physical properties 
of galaxies from multi-wavelength observations and, con- 
versely, to simulate galaxy UV-to-IR fluxes and colours from 
the predictions of theoretical models, both semi-analytic 
and numerical simulations. Semi-analytic models (SAMs, 
Bower et al. 2006; Croton ct al. 2006; Dc Lucia et al. 
120061 ; iMonaco et "all 120071 ; ISomerville et all l20Qgft simulate 
the formation and evolution of galaxies, moving from a 
statistical description of the evolution of the Dark Mat- 
ter structure ("merger trees" either analytically determined 
or extracted from the predictions of numerical simulations) 



and using a well defined set of simple, but physically mo- 
tivated, analytic "recipes" to describe the evolution and 
interplay of the various processes acting on the baryonic 
component (i.e. gas cooling and infall, star formation, stel- 
lar and AGN feedback). SAMs provide detailed informa- 
tion on the physical properties of model galaxies (in terms 
of star formation and AGN activity, stellar mass and gas 
mass content, stellar and gas metallicity) and some infor- 
mation on the relative geometry of the system (the sizes of 
the disc and b ulge components). Stellar popula t ion synthe- 



sis codes (e.g., Fioc fc Rocca-Volmerangg 199 



1 19981 ; iLeitherer et all II 1 19991 ; iBruzual 



.Silva et al.l 
2003h . 



Chariot! I2UU3I ). cou 
pled with the star formation histories predicted by the SAM, 
may be used to produce synthetic spectra and photometry 
for the unattenuated starlight. However, most SAMs do not 
follow the evolution of the relationship between dust mass 
and gas and metal content, nor the relative geometry be- 
tween stars and dust, nor the detailed dust properties (grain 
size distributio n and composition). In the past years, several 
groups (see e.g Deyriendt fc Guiderdoni 2000] ; Hatton et al 



2003 



2007 



Baugh et al.l 120051 ; ISilva et al.l 120051 ; iFontanot et al. 
Lacev et alj l2008f ) have presented SAM predictions 



with refined treatments for dust attenuation and emission 
based on different assumptions for dust grain properties and 
the relative star/gas/dust geometry. 

The various spectrophotometric codes presented in the 
literature show reasonable agreement on the intrinsic (pure- 
stellar) galactic SEDs, but they differ in the treatment of 
the highly complex effects of dust. Assuming that all energy 
absorbed by the dust at UV and optical wavelengths is re- 
radiated in the IR, it is possible to predict the expected IR 
fluxes in the region of interest. The approaches that have 
been used to date can be divided into two broad categories: 
those that use codes explicitly solving the equations of ra- 
diative transfer (RT) in a dusty medium, and those that 
make use of analytically motivated and/or empirically cal- 
ibrated IR template libraries. Clearly the two approaches 
differ in the complexity and (therefore) in the computa- 
tional effort involved in the treatment of dust effects. IR 
templates are usually b uilt up by combin i ng the results 
of analytic calculations |Desert et al.1 Il990l ; iDale fc Heloul 
|2002| ; lLagache et aD |20o4T and/or available observational 
constraints (|Charv fc Elbazl 1200 ll ; iRieke et all 120091 ). They 
are then calibrated by means of comparison with a well de- 
fined set of low-z objects. It is therefore unclear to which 
extent the same templates provide a good description of 
the SED properties a t higher redshifts: several authors (see 
e.g IRieke et al l l2009h discuss possible limitations in using 
these templates outside the redshift range on which they 
are calibrated. RT-solver methods of course provide more 
detailed predictions and a more correct dependence on the 
physical properties and geometrical configurations of dust 
in galaxies. However, the use of a full RT-solver coupled to 
a SAM substantially reduces the efficiency and flexibility of 
the semi-analytic appr oach (becoming the b ottleneck of the 
computation, see e.g. IFontanot et al. [2007). Furthermore, 
this approach requires, and is sensitive to, a large number 
of additional parameters (specifying the details of the dust 
model) which must still be chosen somewhat empirically and 
which carry along numerous assumptions (such as constancy 
with redshift and/or environment). Given these drawbacks, 
and bearing in mind the large number of uncertainties and 



Dust modeling in SAMs II 3 



approximations already inherent in the SAM modeling, it is 
worth asking whether the difference between the results ob- 
tained with a full RT treatment versus IR templates justifies 
the use of the former tool. 

This paper is the second in a series aimed at better 
understanding and characterising the different approaches 
for determ ining dusty SEDs in the SAM framework. In a 
first paper l|Fontanot et aljfeoosl . hereafter Paper I) we stud- 
ied the effect of different prescriptions for dust attenuation 
at ultra-violet and optical wavelengths. In this paper we 
will instead focus on the re-emission of the absorbed en- 
ergy in the infrared region (from 2/im to the sub-mm). In 
the light of the previous discussion, we set out here to ad- 
dress a number of issues. First we briefly recap the study 
that we carried out in Paper I, in order to understand how 
well the amount of light absorbed by dust is predicted by 
analytic recipes for dust attenuation compared with a RT- 
solver. Our main goal in then to compare, on a statistical 
basis, the SED shapes contained in various IR template li- 
braries from the literature with the predictions of a full RT 
solver coupled with a SAM and dust model. Both aspects 
of SED modeling will be of fundamental importance, given 
the growing interest triggered by infrared observations and 
the future availability of large datasets, thanks to facilities 
like the Herschel satellite. As a final goal for this project, 
we identify critical spectral regions, where the choice be- 
tween RT-solver and IR templates may bias the comparison 
of model predictions with available observational data. To 
achieve these goals, we make use of the GRASIL RT code cou- 
pled with the MORGANA semi-analytic model. We derive a set 
of IR SED templates from the coupled MORG AN A+ GRASIL 
outputs, binned as a function of bolometric luminosity, and 
compare these with a range of templates from the literature. 
Finally, we compare the IR luminosity functions predicted 
by the full MORGANA+GRASIL runs with those that we obtain 
if we substitute the template approach. It is worth stress- 
ing that GRASIL is of course not the only RT-solver available 
in the literature. Several other groups have developed al- 
ternative approaches and techniques (se e e.g IPopescu et al.l 
120001 ; iDopita et ail 120051 ; IJonssonl l200rj ) for the prediction 
of panchromatic SEDs, based on different assumptions for 
the properties of star forming regions, dust composition and 
distribution as well as for the star/gas/dust geometry. A 
detailed comparison of the performances of these codes, ap- 
plied to the same sample of star formation histories, would 
be extremely interesting but is clearly beyond the scope of 
the present paper. Nonetheless, in the following we will refer 
to the MORGANA+GRASIL predictions as being representative 
for the performance of a generic SAM+RT solver code, and 
leave further comparison between the different approaches 
to future work. 



This paper is organised as follows. In sec. 12.11 we de- 
scribe the main features of the MORGANA and GRASIL mod- 
els and we review the procedure used in Paper I to define 
the synthetic RT-based SED libraries that we will use in 
the comparison. We introduce the IR template libraries in 
sec. 12.31 and we summarise the main results of Paper I in 
sec. 12.41 In sec. I3.2l we present the statistical comparison be- 
tween the properties of the RT-based SEDs and the IR tem- 
plates, while in sec. 13.31 we discuss the effect of the different 
prescriptions on the statistical properties of model galaxies 



through the analysis of IR luminosity functions. Finally in 
section [4] we present our conclusions. 



2 GALAXY MODELS 

2.1 Semi-analytic model: MORGANA 

We make use of the semi-ana lytic model MORGANA , devel- 
oped by (|Monaco et all 120071) and iFontanot et all l|2007h . 
We refer the reader to these papers and Paper I for more 
details on the galaxy formation model, and here we just re- 
call its main features: the model implements a sophisticated 
treatment of mass and energy flows between the different gas 
phases (cold, hot and stars) and galactic components (bulge, 
disc and halo), as well as a new treatm ent for gas cool- 
ing and infall (following IViola et al.l [200^ 1. It also includes 
both a multi-p hase descriptio n of star formation and feed- 
back (following lMonacoll2004r ) and a self consistent descrip - 
tion of AGN activity and feedback (IFontanot et al.l [2006) . 
For consistency with Paper I, we use the same star for- 
mation historie s (SFHs) and SED s ample Jj], extracted from 
the MORGANA fMonaco et al.|[2007T ) realization presented in 
IFontanot et al.l ( 20071 ). In this study we are mainly inter- 
ested in the IR emission due to the coupling of stellar ac- 
tivity with the dusty interstellar medium, and we do not 
explicitly include the contribution of the central AGN to 
the predic ted SED s . This particular MORGANA realization 
assumes a lSalpete 3 j 19551 ) Initial Mass Function with mass 
range from 0.1 to 100 M©. 

Every model galaxy is represented assuming a compos- 
ite geometry including both a spheroid and a disc compo- 
nent. Dis c expo nential profiles are computed following the 
I Mo et all (|l998l ) formalism: the spin parameter of the DM 
halo is randomly extracted from a well defined distribution 
and the angular momentum is conserved. Bulge sizes are 
computed assu ming that the kin etic energy is conserved in 
merger events |Cole et al.ll200ol ). The presence of a bulge 
is taken into account when disc sizes are computed. MOR- 
GANA provides predictions for the star formation history, 
metal enrichment, and mass assembly of each component 
separately. Thi s information is t hen interfaced with the RT- 
solver GRASIL (ISilva et al.lll998T ) in order to predict the re- 
sulting SED from the UV to the Radio. The MORGANA real- 
ization we consider is able to reproduce the local and 2 = 1 
stellar mass function, the cosmic star formation history, the 
evolution of the stellar mass density, the slope and normal- 
isation of the Tully-Fisher relation for spiral discs, the red- 
shift distribution and luminosity function evolution for K- 
band selected samples, and, more interestingly, the number 
counts of 850 micron selected sources. Despite these suc- 
cesses, the agreement between model predictions and ob- 
servations of the apparent "dow nsized" trend of galaxy for- 
mation is still under debate (see IFontanot et"aL I l2007l , 120091 
for a complete discussion about "downsizing" trends and 
SAMs). It is well established that this model overpredicts the 



1 In Paper I we refer to these ensembles as ML libraries. For 
the sake of simplicity and to avoid further confusion with the 
definition of the IR template libraries, in the following we will 
simply refer to the SFHs we extract from MORGANA as samples 
(i.e. low-z and high-z samples). 
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Figure 1. Left panels: absorbed starlight LiRas a function of stellar mass for the objects in our samples. Green dots refer to the GRASIL 
results, while red circles with errorbars represent the mean predictions of the 4 prescriptions presented in Paper I, as indicated in the 
labels. Right panels: residuals in the relations (symbols and labels as in the left panels). 



number of faint, low-mass galaxies at z < 2 (|Fontana et al.l 
l200fj ; iFontanot et al.ll2009h and the spa ce density of bright 
galaxies at z < 1 ( Monaco et al] 12006). Moreover, it does 
not reproduce the observed levels of star formation activity 
as a function of st ellar and halo mass in the local universe 
(|Kimm et al]|2009h . 

Many relevant properties of galaxies show redshift evo- 
lution, including star formation rates, sizes, metal and dust 
content. In order to test the effects of the predicted changes 
in these physical quantities, we follow the same approach as 
in Paper I, i.e. we draw two different samples of SFHs from 
the MORGANA predictions: a low-z (0.0 < z < 0.2) and a 
high-z sample (2 < z < 3). In Paper I, model galaxies in 
each sample are also split into disc-dominated and bulge- 
dominated subsample^j to test the effect of the composite 
geometry. We check that the shape of the IR SED does not 
show any strong dependence on the details of the geomet- 
rical configuration and our results are not affected by the 
splitting. Unless otherwise explicitly stated, in the following 
we will lift this subdivision. 



2.2 RT solver: GRASIL 

For each object in our samples we interface MORG ANA predic- 
tions with the spectro-photometric code GRASIL l|Silva et al.l 
1 1994 fo r sub s equent improvement we refer the reader 



to 



Silval Il999[: iGranato et all l2000l : iBressan et ail |2002t 



iPanuzzo et al 



2003l : IVega et al , to compute the corre- 



sponding SEDs from the UV to the Radio. GRASIL solves the 



equation of RT, taking into account a state-of-the-art treat- 
ment of dust effects. Stars and dust are distributed in a bulge 
(King profile) + disc (radial and vertical exponential pro- 
files) axisymmetric geometry. The clumping of both (young) 
stars and dust through a two-phase interstellar medium with 
dense giant star-forming molecular clouds embedded in a dif- 
fuse ("cirrus") phase are considered. The stars are assumed 
to be born within the optically thick MCs and to gradually 
escape from them on a time-scale t csc , giving rise to age- 
wavelength-) dependent extinction (i.e the youngest and 
most luminous stars suffer larger extinction than older ones). 
The dust composition consists of graphite and silicate grains 
(with a distribution of grain sizes), and Polycyclic Aromatic 
Hydrocarbons (PAH) molecules. At each point within the 
galaxy and for each grain type the appropriate tempera- 
ture T is computed (either the equilibrium T for big grains 
or a probability distribution for small grains and PAH^f]). 
The radiative transfer of starlight through dust is computed 
along the required line of sight yielding the e merging SED. 
The simple stellar population (SSP) library l|Bressan et al.l 
1 19981 . 12003 ) includes the effect of the dusty envelopes around 
AGB stars, and the radio emission from synchrotron radia- 
tion and from ionised gas in HII regions. Although GRASIL 
computes the resulting SEDs along different lines of sight, 
in the following we use only angle- averaged SEDs. 

Most of the physical information (such as stellar mass, 
star formation history, cold gas mass, gas metallicity) and 
parameters (scale radii for stars and dust in the disc and 
bulge components) needed by GRASIL are provided directly 
by MORGANA. We fixed the parameters used by GRASIL and 



2 according to their bulge-to-total ratio, with a threshold value 
of 0.6 



3 The detailed P AH emission spectrum has bee n updated in 
IVega et al] ||2005| ) based on the lLi fc Draind ll200ll ) model. 
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not provided by MORGANA as in iFontanot et alj (120071 ) and 
Paper I. (i) The disc scale heights for stars and dust are set to 
0.1 times the corresponding scale radii, (ii) We set the escape 
time-scale of young stars for the parent MCs to t csc = 10 7 
yr. This is a good compromise between the value needed 
to explain the SE P of spirals (~ a few Myr) and starbursts 
{Silva et al.lll998l . ~ a few 10 Myr, see) and it is similar to the 
est imated destruct i on tim e scale of MCs by massive stars. 
In IFontanot et al.l {2007) this value was found to produce 
good agreement with the K-band luminosity functions and 
the 850 ^tm counts, (iii) The total gas mass is split between 
the dense and diffuse phases, assuming that 50% of the gas 
is in the star forming molecular clouds (/mc = 0.5). The 
results are not very sensitive to this choice, (iv) The mass of 
dust is obtained by the gas mass and the dust-to-gas mass 
ratio Jdust which is set to evolve linearly with the metallicity 
given by the galaxy model, <Wst = 0.45 Z. (v) The optical 
depth of MCs depends on the ratio tmc oc 8 dust Mmc/>"mc; 
we set the mass and radius of MCs to typical values for the 
Milky Way, M M c = 10 6 M Q and r M c = 16 pc. (iv) The dust 
grain size distribution and composition is chosen to match 
the mean Milky Way extinction curve. 

Many authors have suggested that dust properties may 
evolve as a result of the evolution of galac tic properties 
and / or differential metal enrichment (s e e e.g. ISchurer et al.l 
2009). In a recent paper ILq~ Faro et all (|2009l ') analysed the 
luminosity function of Lyman Break galaxies at z > 3 in the 
MORGANA framework, assuming the same GRASIL choices for 
the dust properties and their relation with gas and metal 
content. Their results show th at the predicted exti nctions 
are larger than the estimates of lBouwens et all l|2007l ). They 
then suggest that different values for t csc and /mc are able 
to reproduce the correct amounts of attenuation. Given the 
great uncertainties in the expected properties of dust as a 
function of redshift, we assume the same relations between 
dust and gas content, as well as the same dust composition 
and grain distribution for both the high-z and low-z sam- 
ple. This assumption is similar in spirit to applying an IR 
template library derived from low-redshift observations at 
higher redshift. Thus, any differences that we find in the 
IR SEDs as a function of redshift are due to the changing 
physical properties of the model galaxies, not to changes in 
the dust properties. Given the strong dependence of the pre- 
dicted SEDs on the details of the SFHs of model galaxies, 
the quantitative results we present in this paper are strictly 
valid only for the MORGANA+GRASIL algorithm, i.e. by in- 
terfacing GRASIL with another SAM, we expect to obtain a 
different sample of synthetic SED libraries. Nonetheless, the 
comparison between the predictions of these particular SED 
libraries and IR templates will provide us fundamental in- 
sight on the dependence of SAM pre dictions on the differen t 
IR modeling. Moreover, as shown in IFontanot et~al] l|2009h . 
the global evolution of bulk galaxy properties (such as stel- 
lar masses and star formation rates) in GRASIL agrees well 
with other SAMs in the literature, suggesting that the re- 
sults presented here would likely be similar for any currently 
available SAM. 



2.3 IR template libraries 

Starting from the pioneering work by iDesert et al. (1990, 
hereafter DSP90), many authors proposed analytic solutions 



to RT equations, and used their results to define template 
IR SEDs as a function of Lir, calibrated by comparison 
with known local prototypes. The original formulation of 
DSP90 was meant to reproduce the dust emission from indi- 
vidual regions within the Milky- Way. They assumed 3 main 
contributors to dust emission at IR wavelengths: polycyclic 
aromatic hydrocarbons (PAHs), very small grains and big 
grains. The former two are composed of graphite and sil- 
icates, with small and big grains probably dominated by 
graphite and silicate respectively. The thermal properties of 
each species are defined by its size distribution and its ther- 
mal state. Big grains are assumed to be in near thermal equi- 
librium. Their emission can be modeled as a modified black- 
body spectrum. On the contrary, small grains and PAHs are 
very likely in a state that is intermediate between thermal 
equilibrium and single photon heating. They are therefore 
subject to temperature fluctuations and their emission spec- 
tra are much broader than a modified black-body spectrum. 

In the original DSP90 paper, this approach was suc- 
cessfully applied to rep roduce the SED of the Milky Way. 
iDevriendt et al.1 (|l999l . hereafter DGS99) expanded this 
framework to define IR template libraries: since the de- 
tailed size distributions are modeled using free parameters, 
it is possible to calibrate them by requiring the model to 
fit a series of observational constraints, such as the extinc- 
tion/attenuation curves, observed IR colours and the IR 
spectra of local galaxies. Once the free parameters are fixed 
it is possible to use this model to predict the IR spectral 
contribution corresponding to each species, embedded in ra- 
diation fields of different intensity. The distribution of dust 
mass over heating intensity is usually assumed to follow 
a power-law. The main advantage of this approach lies in 
the relatively small number of parameters required by the 
model (with respect to a RT solver). On the other hand, it 
relates a single SED to a class of model galaxies (defined 
by their infrared luminosity -Lir), irrespective of their mor- 
phologies or other pr operties. Alternative implementations 
have been proposed bvlDale fc Heloul (|2002l , hereafter DH02) 
and lLagache et al.l l|2004l . hereafter L04) . These implemen- 
tations differ from the original DSP90 work, mainly in the 
modeling of the radiation field (DGS99, DH02), and in the 
relative contribution and shape of the IR emission of the 
different species (L04, DH02). The DGS99 template library 
consists of nine SEDs, corresponding to values spanning 
the range 10 6 L Q < £ir< 10 14 L©; L04 consists of five SEDs 
in the range 10 9 Lq < -Lir< lO 13 !/©; DH02 contains sixty- 
four SEDs in the range 2 x 1O 8 L < L m < 2 x 10 14 L Q . 

An alternative approach has been proposed by 
ICharv fc Elba3 (|200ll . hereafter CE01). They start from four 
reference SEDs generated using GRASIL and constrained to 
reproduce the observed SED of Arp220, NGC6090, M82 an d 
M51 from the far-UV to the sub- mm j|Silva et all ll998l V 
These four galaxies are considered representative of four 
galactic populations or luminosity classes: ULIRGs, LIRGs, 
starburst and normal galaxies respectively. They calibrate 
the 3(j,m < A < 18/im region of the four model spec- 
tra with ISOCAM observations, then they split each SED 
into a mid-IR and a far-IR section and interpolate between 
the reference spectra to generate a set of sample template 
SEDs at intermediate luminosities. Additional templates 
from iDale et all j|200ll ) are used to widen the range of spec- 
tral shapes. The final templates are then chosen by requiring 
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Figure 2. IR templates from the literature (red solid lines) compared to the GRASIL predicted mean SEDs for the CE01, DH02 and 
L04 libraries and four representative values of Lir. In each panel the dashed and dotted lines refer to the low-z and high-z sample 
respectively. It is worth noting that the apparently different shapes of the mean SEDs in the various panels are due to the different 
LiRbinning considered in the original libraries (see text for more details). 



the two sections to account for a variety of observational con- 
straints and merging together the best fit solutions in the 
two regions. The CE01 library consists of 105 templates, 
covering the range 2.7 x 10 8 L Q < L m < 3.5 x 10 13 L Q . 

Finally, we also consider the recent work of lRieke et al.l 

(|2009l . hereafter R09). They assemble detailed SEDs for 
eleven LIRGs and ULIRGs from observational datasets (in- 
cluding published ISO, IRAS and NICMOS data as well as 
previously unpublished IRAC, MIPS and IRS observations) 
and use them to build representative average IR templates. 
In regions where no hom ogeneous spectral covera ge is avail- 
able, they use GALAXEV (|Bruzual fc Cha riot 2003) synthetic 
spectra, calibrated by requiring them to fit the available 
photometry. To model the far infrared SEDs, they assume 
a single black body with wavelength-dependent emissivity, 
while in the Radio they use a single power law wavelength 
dependence, as suggested by observed data. The R09 li- 



brary includes fourteen SEDs covering the 5.6 x 10 9 L Q < 
Lir< 1O 13 L range. 

2.4 Analytic Models of Dust Attenuation 

In Paper I, we considered the effect of dust at optical to 
near-infrared wavelengths as predicted by GRASIL, by com- 
paring the synthetic SED to the intrinsic starlight emission 
of each model galaxy. We compute mean dust attenuations 
and optical-to-near-infrared colours in the two samples and 
we then compare them to the predictions of simpler ana- 
lytic prescriptions, commonly used in the SAM framework. 
These prescriptions are based on (i) an analytic attenua- 
tion law, normalised to (ii) an assumed value for the face-on 
optical depth at a reference wavelength (usually in the V- 
band, ry), eventually computed as a function of the phys- 
ical properties of the galaxy itself; (iii) an inclination cor- 
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Figure 3. R09 IR templates (red solid lines) compared with the GRASIL predicted mean SEDs. In each panel the LiR,reference interval 
is indicated, while dashed and dotted lines refer to the low-z and high-z sample respectively. 



rection, performed by means of simple analytic models rep- 
resenting simplified geometries (such as the slab model or 
the oblate ellipsoid). In Paper I we tested different choices 
and combinations of these three elements. In parti cular, 
we considered : (i) th e Milky Way jMathis et alj 1 1983ft. the 
ICalzetti et al.l (|200dh and the age dependent Chariot fc Falll 



(2000) attenuation laws, p l us th e composit43 attenuation 



( zIUUU ) attenuation laws, pl us th i 
law of |Pe Lucia fc Blaizotj (l2007h 



(ii) the fitting fo r mulae 
for ry proposed bv | Guiderdoni fc Rocca-Volmerangel (Il987l ) 
and |Pe Lucia fc Blaizotl (|2007l ) ; (iii) slab and oblate ellip- 
soid geometries. Our results show that GRASIL RT calcu- 
lations can be reproduced with reasonable agreement even 
by a model as simple as a slab geometry combined with 



i.e the Milky Way extinc tion curve descr i bes t he effect of the 
diffuse "cirrus" dust and the lCharlot fc Falll d200Cl> power-law the 
attenuation law experienced by young stars in the dense birth 
clouds. 



the composite age-dependent attenuation curve suggested 
by |De Lucia fc Blaizotj [20071 ). if an accurate estimate of 
the intrinsic tv is used to normalise the attenuation law. 
Using these synthetic SED samples, in Paper I we tested 
different analytic expressions, relating tv with the physi- 
cal properties of the model galaxy (such as bolometric lu- 
minosity, stellar mass, metallicity, cold gas fraction and/or 
surface density). We demonstrated that the intrinsic value 
for tv, as predicted by GRASIL, correlates strongly with 
the gas metallicity, the cold gas mass and the scale ra- 
dius of the model galaxies and we provided simple an- 
alytic fitting formulae. We then compared our formulae 
with analogous prescripti ons found in the literature, and 
in par ticula r with the the iGuiderdoni fc Rocca-Volmerangel 
l|l987h and |Pe Lucia fc Blaizotl (|2007f ) fitting formulae, by 
defining four different analytic prescriptions for dust at- 
tenuation and testing them against GRASIL predictions. 
The four prescriptions share the same choice for at- 
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tenuation law and geometry (i.e. a |Pe Lucia fc Blaizotj 
composite attenuation law combined with a slab 
model), but assume a different normalisation (i.e. tv). A 
first prescription uses the intrinsic tv value as computed 
from the comparison of attenuated and unextinguished 
GRASIL synthe tic SEDs (TAU-GS). In the remaining w e as- 
sume both the | Guiderdom fc Rocca-Volmerangel (|l987T ) and 
|Pe Lucia fc Blaizotj ( 2007 ) ry fitting formulae plus our re- 
sults: we then define TAU-GRV, TAU-DLB and TAU-FIT 
prescriptions respectively. We showed that the agreement 
among the different prescriptions is satisfactory for the low- 
z sample, but the discrepancies are significant in the high-z 
sample. Our proposed fitting formulae provide the best com- 
promise in reproducing the GRASIL RT-predictions in both 
samples at the same time. 



3 RESULTS 

3.1 Total Absorbed Starlight 

In this paper we extend the analysis of Paper I by con- 
sidering the energy re-emitted by dust at longer wave- 
length. Therefore, as a first critical test we compare the 
total amount of absorbed starlight from the UV to the near- 
infrared (Lir) both in GRASIL and in the four analytic ex- 
tinction prescriptions. As described in the previous section, 
in defining our prescriptions we assume a fixed combination 
for the geometry and the attenuation law. It is worth stress- 
ing that alternative analytic dust models presented in the lit- 
erature do not necessary assume fixed attenua t ion la ws: for 
example in iGuiderdoni fc Rocca-Volmerangel l)l987h a de- 
pendence of the shape of the attenuation law on metallicity 
is explicitly considered. We compute the intrinsic value of 



Ljfby comparing the extinguished F£ ba and unextinguished 



r GS 

Mr 



SEDs: 



/ r-iunabs pa6s\ t\ 



(i) 



Starting from F™ na s and applying the four prescrip- 
tions for dust attenuation, we estimate the corresponding 
LiRvalues. We compare them with L^ffin fig.[U in each case 
we see reasonable agreement between the predicted total 
amounts of absorbed starlight and GRASIL results. This en- 
couraging result is not obvious, since our prescriptions as- 
sumed a fixed shape for the extinction law (even if we con- 
sider an age-dependent component). This implies a fixed ra- 
tio between the attenuation at different wavelengths. On the 
other hand, we showed in Paper I (fig. 3) that the scatter 
in GRASIL predicted attenuation laws is significant, and al- 
though the Milky Way attenuation law is a good estimate 
for the mean shape, single objects may deviate considerably. 
The agreement is especially good for the disc-dominated 
objects, where we see a tight relation between the stellar 
mass and £§1 • The intrinsic scatter increases in the bulge- 
dominated objects and the analytic prescription are able to 
reproduce the general trend, even if the residuals of the re- 
lation (right panels) show a 0.5dex scatter. This result sug- 
gests that the deviation of synthetic SEDs from a Milky Way 
attenuation law is larger in bulge-dominated objects. Under 
the hypothesis that all the absorbed energy is re-emitted at 
infrared wavelengths, and assuming a set of IR templates as 
described above, it is then possible to estimate the corre- 
sponding infrared SEDs. 



Dust modeling in SAMs II 9 



3.2 Comparing IR SEDs 

In order to compare the predictions of the full GRASIL RT 
calculation to the IR template approach we split the model 
galaxies in the low-z and high-z samples into IR luminosity 
classes, according to their L§f . For each IR template library, 
we define luminosity classes by sorting each model galaxy 
into the nearest LiRbin. We renormalise each GRASIL SED 
in our sample to the ratio between Lf]f and the Lmof the bin. 
We then compute the statistical properties of the synthetic 
SEDs for each class: we create a mean SED, by considering 
the mean flux at each wavelength, and we define a scatter 
in the flux distribution as a function of wavelength. For the 
CE01 templates we consider both the contribution of direct 
starlight and dust emission, whereas for the other libraries 
we consider only the contribution of the dust emission. We 
consider the low-z and high-z samples separately. 

We compare the mean GRASIL SEDs to the correspond- 
ing IR templates, for four representative values of LiRand 
three libraries in fig. [2] The agreement between the mean 
SEDs constructed from our samples and the IR templates 
is satisfactory in most cases, but there are some significant 
discrepancies. In order to perform the same comparison over 
a wider range of Lir, we consider the full dynamical range 
probed by the R09 library in fig. [3] The most prominent 
difference between the templates and the predictions of the 
GRASIL RT-solver concerns the position of the peak of the 
thermal dust emission: the peak is shifted to longer wave- 
lengths in the GRASIL predictions relative to the IR tem- 
plates. It is particularly evident in fig. [3] that the discrep- 
ancy depends also on Lir. In order to get a better insight 
into this problem, we compare the mean GRASIL SEDs for 
different values of the total amount of energy emitted in 
the IR in fig. [3] (left panel) . Lmhas a strong effect on the 
overall normalisation of the mean SEDs, but they show a 
quite similar shape. In particular, the position of the IR 
peak does not depend on Lir over the whole range studied, 
while it is evident from fig. El that it is changing significantly 
in the templates. Another notable discrepancy is present in 
the Mid-Infrared (8/im< A < 20/im), where the highly un- 
certain contribution of PAH emission is dominant. Finally, 
it is also worth noting the strong discrepancy seen in most 
cases at the shortest wavelengths (A < 8/im). This is the 
region where direct starlight emission from old stars is still 
significant in the galactic SEDs and comparable to the dust 
thermal emission: therefore we interpret this result as due 
to the difficulties in disentangling the relative contribution 
of these two components. 

The discrepancy between the IR templates and the 
mean SEDs for the low-z sample is particularly interest- 
ing, since the former have been calibrated to the proper- 
ties of local galaxies, and since GRASIL is able to repro- 
duce thei r SEDs, under a reasonable choice for dust pa- 
rameters (ISilva et al.lll998T ). In particular, the position de- 
pends mainly on the temperature distribution of the emit- 
ting grains, and it is likely connected to the assumed model- 
ing of dust abundance and distribution. This discrepancy is 
due to the combination of two separate and degenerate ef- 
fects. First of all, we assume the same dust properties for all 
galaxies in our samples. This represents a common choice 
in the semi-analytic framework, and it provides a reason- 
able description for large model galaxy samples. Second, our 



synthetic SEDs depend on the MORGANA predicted physical 
conditions o f model galaxies; we already know from pre- 
vious work (|Fontanot et all 120071 ; iKimm et all 120091 ') that 
this particular MORGANA realization does not reproduce all 
properties of local galaxies, i.e. their levels of star formation 
and their relative gas and stellar content. Most notably, in 
the original formulation of the model we do not attempt to 
calibrate our prediction to reproduce the IR properties of 
local galaxies. Therefore it is perfectly plausible to obtain 
mean SEDs which differ from the local sample. Nonetheless, 
the proposed comparison is instructive regarding the differ- 
ent results obtained by the two approaches for the generic 
galaxy sample extracted from a cosmological galaxy forma- 
tion model. 

The most interesting insight from the analysis of fig. [5] 
and [3] is the different shapes of the mean SEDs drawn from 
the high-z and low-z sample. In most cases the mean SEDs 
corresponding to the high-z sample show on average a better 
agreement with the corresponding template. This result is 
somewhat surprising, given the fact that IR templates are 
calibrated using low-z observations. In order to better under- 
stand this result, we construct the mean SEDs over a finer 
binning in redshift, while retaining the same LiRbinning as 
for the R09 templates. We compare the evolution of the new 
mean SED sample with redshift in fig. [4] (right panel), for the 
representative luminosity bin centred at Log(Lm) = 10.5. 
Different luminosity classes show similar results. It is evi- 
dent from this plot that the mean SED is experiencing a 
strong redshift evolution, in particular in the spectral re- 
gions corresponding to the peak of thermal IR emission and 
PAH emission. This pattern can be ascribed to the redshift 
evolution of the physical properties (i.e. gas content, metal 
enrichment and size evolution) of the model galaxies in the 
MORGANA realization and represents a warning against us- 
ing the same template library over a wide redshift range (as 
already noticed by DN02 and R09) . It is also worth mention- 
ing that the same dependence of the mean SED on Lir at 
fixed redshift (fig.[H left panel) seen at z = holds at higher 
redshifts. In order to get better insight into the main driver 
of the evolution of the SED shape, we perform a fitting pro- 
cedure similar to Paper I (see their sec. 4). For each synthetic 
SED in our library, we modeled the wavelength correspond- 
ing to the peak of IR emission as a power-law function of 
the physical properties of the corresponding model galaxy. 
We consider Lir, the stellar and cold gas masses, the star 
formation rate, metallicity and galactic radius and we define 
a set of general relations involving independent quantities. 
For each combination, we determine the best-fit parameters 
and evaluate the goodness of the fit through a \ 2 procedure. 
Our results show that the strongest correlation is found for 
the total mass surface density of the galaxy. This is not com- 
pletely unexpected since the relative spatial distribution of 
stars and gas is fundamental to predicting the temperature 
of dust grains. The scatter in the correlation is large, but we 
find that adding additional degrees of freedom does not re- 
duce it considerably. We conclude that other physical quan- 
tities, like star formation rate and metallicity, still play a 
non-negligible role in determining the resulting shape of the 
SED. 

We then consider the distribution of fluxes as a func- 
tion of wavelength in fig. [5] It is evident from the plot that 
the scatter within each LiRclass may be important, espe- 
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Figure 5. Left Panels: The 1-cr distribution of GRASIL predicted SEDs in the same LiRreference intervals and for the same libraries as 
in fig[2] Flight Panels: The 1-cr distribution of GRASIL predicted SEDs for the same LjRreference intervals as defined in R09. The central 
Ltr, value for each interval is indicated in each panel. In all panels the yellow shaded area and horizontal texture refer to the low-z and 
high-z sample respectively. 



cially at the lowest values of total absorbed starlight. Both 
the scatter variations between different libraries at a given 
LiR(left panels) and along the same library (right panels) 
are mainly due to the inhomogeneous distribution of galax- 
ies in the LiR,bins. As expected, larger binnings correspond 
to larger scatters around the mean SEDs. At the same time, 
the variance of the flux as a function of wavelength does 
not show any significant dependence on redshift, for a given 
Lmbin. This result implies that using a single SED from 
a template library to describe all galaxies belonging to a 
given LiRclass, may introduce a bias between the compari- 
son between model predictions and observations. Our anal- 
ysis shows that it is possible to use mean SEDs as represen- 
tative for all SEDs belonging to a luminosity class only in a 
statistical sense (and this proxy improves with finer binning 
in I/ir), but not on a object-by-object basis. It is also worth 
noting that the different binning has a non- negligible effect 
on the shape and normalisation of the corresponding mean 
SED, due to the different sample definition. Although the IR 
templates share a similar trend, this is another hint that a 
finer binning helps reduce the differences between templates 
and GRASIL RT predictions. 
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3.3 Effect on the IR Luminosity functions 

We then compare the individual monochromatic luminosi- 
ties for each model galaxy as predicted by the GRASIL com- 
putation and by the IR templates. We select four wave- 
lengths widely used in the literature to track different galaxy 
populations and constrain the properties of dust and PAHs 
at different temperatures: namely the 8fim, the 24^m, the 
60jj,m and the 850/im restframe luminosities. For each model 



Figure 6. Intrinsic IR luminosities at 8, 24, 60, 850 fira compared 
with the corresponding predictions for the R09 templates. 

galaxy we estimate the expected fluxes from the IR template 
library as follows. We use the intrinsic GRASlL-predicted 
L^to select the closest template for each IR library. We 
then rescale the template to the ratio between its LiRand 
LrR , and we use it to estimate the monochromatic luminosi- 
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ties in the four bands. For the CE01 templates we simply 
use the template, while for the other libraries we sum up the 
chosen template with the extinguished starlight spectrum 
(as predicted by GRASIl). We then statistically compare the 
predictions derived using the template libraries with the in- 
trinsic results of GRASIL computations by constructing the 
corresponding luminosity functions. In this way we are able 
to quantify the impact of the different choice for IR tem- 
plates on the statistical properties of galaxy samples, as pre- 
dicted by SAMs. It is worth stressing here that we do not 
require our SAM to reproduce the properties of the local 
Universe. Instead our aim is to check the impact of different 
prescriptions for dust emission on the predictions of SAMs. 

In fig. [6] we show a comparison between the monochro- 
matic luminosities obtained using the R09 templates with 
the intrinsic GRASIL predictions (similar results hold for the 
other templates), while in fig. [7| the resulting luminosity 
function for all considered templates (right panel). The in- 
trinsic luminosity functions for our low-z and high-z samples 
are represented by the solid thick line, while other lines cor- 
respond to the predictions obtained using the four IR tem- 
plate libraries. We find similar results with respect to the 
analogous analysis at optical and near-infrared wavelengths 
presented in Paper I (their fig. 12): both the shape and nor- 
malisation of the luminosity functions are correctly recov- 
ered in most cases. These results clearly show that the use 
of a statistical estimator like the luminosity function, defined 
over a large galaxy sample and including the contribution of 
various luminosity classes, is less sensitive to uncertainties 
in the prediction of individual galaxy properties. Nonethe- 
less, some important discrepancies are seen. In particular, 
in the sub-mm region, GRASIL fluxes are systematically un- 
derpredicted by the templates: this is mainly due to the dis- 
crepancy in the position of the ~ 100/^m peak in the SEDs 
(see fig. [5} , which affects the shape of the SEDs at longer 
wavelengths. In the low-z sample, the flux at 850/^m is un- 
derpredicted only for the brightest sources, while at fainter 
fluxes the agreement is reasonable. However, in the high-z 
sample the templates are not able to reproduce the bright- 
est sources at 8/wn, while the agreement at 24/^m and 60/im 
is again satisfactory. It is possible to understand this pecu- 
liar behaviour by looking at the shape of the brightest mean 
SEDs compared with the corresponding templates in fig. [2] 
In order to assess the impact of these effects in terms of the 
comparison between model predictions and data, we show 
in fig. [7] the available observational constraints for the cor- 
responding luminosity functions. We warn the reader that 
we do not explicitly tune our model to reproduce the cor- 
responding samples; moreover, real data samples do not al- 
ways match the same redshift interval we adopt to define 
our SED libraries. Nonetheless, we include them to guide 
the eye and to provide an hint of the effect of different IR 
modeling in the prediction of the luminosity function. As 
expected, the most interesting results are seen for the 8 /im 
and 850 /im bands. In particular, the 8 \im z ~ 2 luminosity 
function is better reproduced by the GRASIL RT predictions, 
while it is underpredicted by IR template approach. On the 
other hand, the low-z 850 \im luminosity function, which is 
clearly overpredicted in the MORGANA+GRASIL framework, 
is recovered by the same galaxy formation model when the 
IR template approach is adopted. 



4 SUMMARY 

This paper is the second in a series aimed at (i) under- 
standing and quantifying the results of detailed calculations 
of radiative transfer in a dusty medium and (ii) compar- 
ing these results with simpler analytic recipes coupled with 
semi-analytic models and spectro-photometric codes. The fi- 
nal goal of this work is to assess the impact of different mod- 
eling for dust attenuation and emission in the framework of 
semi-analytic models of galaxy formation and evolution. In 
the first paper of the series (Paper I) we studied the effects 
of dust attenuation at optical and UV wavelengths, while in 
this paper we focus on the dust emission in the IR region. 
We use the same star formation history libraries based on 
the MORGANA semi-analytic model and defined in Paper I: 
they include a low-z (0 < z < 0.2) and a high-z (2 < z < 3) 
sample, with the aim of assessing the effect of the evolution 
of galaxy physical properties on the corresponding SEDs. 
We couple the star formation histories and galaxy proper- 
ties to the RT-solver GRASIL, to obtain an estimate of the 
synthetic SEDs from the UV to the Radio. 

As a first check we consider the intrinsic absorbed 
starlight Liaaccording to the GRASIL results and we compare 
it to the predictions of the four prescriptions for dust atten- 
uation at optical to near-infrared wavelengths discussed in 
Paper I. As already shown in Paper I, the scatter in the 
GRASIL predicted attenuation laws is not negligible: there- 
fore we expect discrepancies on an object-by-object basis. 
However, our results show that all prescriptions provide a 
reasonable description of the mean intrinsic Liafor the whole 
population, within the errors shown in fig. [1] We conclude 
that such simplified analytic approaches can be used to ob- 
tain statistically valid estimates for Lir. 

We then consider four different IR template libraries, 
commonly coupled to spectro-photometric models to pro- 
vide predictions f or dust emission in the IR. They include 
both theoretical (iDale fc Heloul 120021: iLagache et all 12004) 
and e mpirical templates ( Chary fc Elbazl 2001 : Rieke et al.l 
2009). In each case, we define LiRclasses similar to those 
used to generate the IR templates and we sort the synthetic 
SEDs for the low-z and high-z samples accordingly. For each 
class, we then construct mean SEDs by computing the mean 
fluxes as a function of wavelength. Our results show that, 
overall, the mean SEDs are in fairly good agreement with 
IR templates over a wide spectral range. Significant discrep- 
ancies arise, however, both at wavelengths longer than the 
~ lOOpim peak, due to the negligible evolution of the peak 
position with respect to the IR templates, and for low IR 
luminosities. These discrepancies are mainly due to the fact 
that the MORGANA+GRASIL code has not been explicitly cal- 
ibrated to reproduce the details of dust and physical proper- 
ties of IR bright galaxies in the local Universe, but it assumes 
mean values for the relevant quantities. 

Moreover, by comparing a low-z and a high-z sample we 
see significant differences in the properties of mean SEDs, 
again growing larger for low values of Lir. In order to gain 
better insight into the origin of these discrepancies, we com- 
pute mean SEDs on a finer binning in the Lir and redshift 
space. We then consider the redshift evolution of the mean 
SED at a fixed Lir (and vice versa). The mean SED shape 
shows strong redshift evolution at fixed Lir, and very little 
evolution as a function of Lir at a given cosmic epoch. We 
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Figure 7. GRASIL predicted IR luminosity functions in the low-z and high-z samples (thick solid lines) compared with the predictions 
using the IR template libraries. Solid, dotted, dashed, dot-dashed and long-short dashed lines refer to GRASI L, CE01 , DH02, L04 and 
R09 t emplates respectivel y. Data points represen t a compi lation of local measu rements fr omlHu anfi et ah] l|200 7l. 8 fim ) , iRodighicro et all 
l|20ld . 8 iim and 24 urn) ISaunders et al] l|l99d . 60 urn), iDunne et a.1.1 J2000I . 850 urn). ISerieant fc Harrisonl <2005t 850 /im) and high 
redshift results from lCaputi et al.l J2007I . 8 urn) . iRodiehiero et al.l (1201(3 . 8 fim and 24 um) iMarleau et al.l J2007I . 24 fim). 



quantify the dependence of the SED shape on the physical 
properties of the underlying model galaxies by fitting the 
position of the peak of IR emission. Our results show that 
the main contributor is the redshift evolution of the spatial 
distribution of stars and gas (total mass surface density) as 
predicted by MORGANA. Both star formation rate and metal- 
licity cause a secondary, but non- negligible contribution to 
the broad scatter in the relation. This result reflects the fact 
that the main factor responsible for the detailed shape of the 
GRASIL predicted IR SED is the mean temperature of dust 
grains, which is a complex function of the distribution of 
dust particles (linked to the cold gas surface density and 
metallicity) and of the radiation field (linked to the stellar 
mass surface density and star formation rate). Finally, we 
also show that the scatter around the mean fluxes is signif- 
icant and it depends on the width of the Ltr binning. For 



a sparse sampling of Lir the mean SEDs give only a poor 
description of the SED variety. 

These results suggest that it may not be valid to use the 
same IR template library for modeling galaxies at different 
redshifts, where their physical properties are likely different 
(as in the case of MORGANA galaxies). Moreover, despite the 
fact that many physical properties of galaxies (such as gas 
and metal content, present and past star formation activ- 
ity, disc and bulge sizes) are predicted by the SAM, when 
solving the equation of RT we still have to make assump- 
tions about the physical state of the dust (composition, grain 
size distribution) and about its distribution relative to stars 
and the interstellar medium. Given the present uncertain- 
ties, we assume that these properties do not change as a 
function of galactic properties nor with cosmic time over 
the redshift range we consider. Our choice reflects a param- 



Dust modeling in SAMs II 13 



eter combination that has been proven a dequate for repro- 
ducing the propertie s of z < 3 galaxies (jSilva et al.l 1 19981 ; 
iFontanot et al.l 120071 ). However, a strong systematic varia- 
tion of these parameters (especially t eac and /mc) at higher 
redsh ifts is expected bot h observationallv l|Maiolino et alj 
|2004 ) and theoretically (|Lo Faro et all I2009T) . We expect 
any variation in th e dust parameters and/or composition 
i|Schurer et al.l l2009) to lead to even larger differences in the 
mean SED shape and normalisation. 

We then consider the IR monochromatic luminosities 
at Sfim, 24/zm, 60/xm and 850/im as predicted by the MOR- 
GANA+GRASIL model and by the four IR template libraries. 
We compare the fluxes on an object-by-object basis and the 
resulting luminosity functions. Although there are large dis- 
crepancies for some individual objects (up to two orders of 
magnitude), the overall statistical agreement is quite good 
in most regimes. There is a systematic discrepancy at 850/im 
for both the low-z and high-z samples, with a systematic un- 
derprediction of the GRASIL RT predicted luminosity func- 
tions. This is mainly related to the discrepant position of the 
~ 100/im peak in the SEDs with respect to the templates. 
The brightest sources in the high-z sample show disagree- 
ment also at 8fim, while the 24/im and 60/xm predictions are 
similar. This is related to the peculiar shape of the high-z 
sample mean SED at this IR luminosity, and, again, implies 
that caution should be used when comparing theoretical IR 
predictions from different methods, even on a statistical ba- 
sis. 

Our results extend those presented in Paper I into the 
IR regime, and our conclusions are quite similar. The level of 
agreement between the predictions of a RT solver approach 
(superior in order to understand the details of dust absorp- 
tion and re-emission in galaxies on a object-by-object basis) 
and more computationally efficient semi-analytic methods 
depends on the spectral regions under analysis. We have 
identified wavebands where the agreement is good (in a sta- 
tistical sense) and others where the discrepancies are sig- 
nificant. In particular, for the MORGANA+GRASIL model the 
agreement between the two approaches is good in the Mid- 
IR, while we found significant discrepancies at 8/xm and 
850/im. 

As a wealth of multi-wavelength observations are be- 
coming available in the IR, we expect to use them to con- 
strain with greater accuracy the redshift evolution of dust 
properties and their relation to the physical properties of 
galaxies, as well as the observed shape of their SEDs over a 
wide range of redshift. It has been recently shown that the 
computational time required for a full RT solver coupled to 
a SAM can be considerably reduced , with no loss of accu- 
racy, with the use of neural networks (jSilva et alj201lf ) , thus 
helping to overcome one of the biggest drawbacks connected 
to the use of this tool in the SAM framework. At the same 
time, improvements in the SAMs are needed, i.e. the model- 
ing of dust generation, evolution and dispersion, in order to 
reduce the number of externally fixed parameters and make 
use of the full predictive power of RT solvers. 
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